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Abstract 


Introduction 


The SOUSSA {^teady, _Osci 1 1 atory , and 
Unsteady ^ubsonic and Supersonic Aerodynamics) 
program is the computational implementation of a 
general potential -flow analysis (by the Green's 
function method) that can generate pressure dis- 
tributions on complete aircraft having arbitrary 
shapes, motions, and deformations. This paper 
presents results of some applications of the 
initial release version of this program to 
several wings in steady and oscillatory motion, 
including flutter. The results are validated by 
comparisons with other calculations and 
experiments. Experiences in using the program 
as well as some recent improvements are 
described. 
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generalized aerodynamic force 
wing root semi chord 
lift-curve slope 


pressure coefficient 
lifting-pressure coefficient, 

Cp - Cp 

reduced frequency, y- 


freestream Mach number 

total number of panels on one- 

quarter of wing surface, e.g., on 

upper right-hand side 

Reynolds number based on average 

chord 

freestream speed 

local chordwise coordinate measured 
from leading edge, fraction of local 
chord 

fuselage longitudinal coordinate 
measured from nose, fraction of 
fuselage length 

spanwise coordinate, fraction of 
semi span 

real part, imaginary part 
angle of attack 

mass ratio, mass of wing divided by 
mass of a volume of air at freestream 
density contained within a cylinder 
ci rcumscribed about wing planform 
circular frequency of wing 
osci 1 1 at ion 

circular frequency of first torsion 
mode 


Subscripts: 


u,Ji upper surface, lower surface 


Application of a generalized Green's func- 
tion method to the full, time-dependent 
potential -flow equation leads to an integral 
equation for the velocity potential at any point 
in the flow, including points on the surface of 
a body or bodies in the flow. The SOUSSA 
(Steady, Oscillatory, and JJnsteady _SiJbsonic and 
S^upersonic Aerodynamics) Pl.l* computer 
program^*^ is a panel-method code which 
implements this integral equation for linearized 
subsonic flow in the complex-frequency domain, 
and within that context is applicable to general 
shapes such as complete aircraft or other bodies 
having aribitrary shapes, motions, and 
deformations. A program with this degree of 
generality has many possible uses, including (1) 
unsteady-state applications, such as flutter, 
gust-response, and active-controls analyses with 
multiple sets of frequencies, mode shapes, and 
Mach numbers, and (2) steady or quasi -steady 
state applications such as aerodynamic analyses 
requiring pressure distributions and aerodynamic 
coefficients, static or dynamic stability analy- 
ses requiring stability derivatives (including 
rate derivatives), static aeroelastic analyses, 
and structural load calculations. The purpose 
of this paper is to present the results of some 
applications of SOUSSA Pl.l to several wings in 
both steady and oscillatory motion, including 
flutter, to validate the results by comparisons 
with other calculations and experiments, and to 
describe experiences in using the program as 
well as some recent improvements made to it. 

The SOUSSA panel method was not formulated 
primarily for application to isolated wings of 
simple shape, nor was it intended to be a com- 
petitor of lifting-surface theory. However, for 
validation purposes the present study focuses on 
pressure-distribution and flutter calculations 
for such simple shapes so that comparisons can 
be made with both steady and unsteady lifting- 
surface calculations as well as with existing 
experimental data. Accordingly, the present 
study i ncl udes ■ aerodynamic calculations for two 
rectangular wings, b a clipped-tip delta 
wing, 6 and two swept wings. Calculations were 
made for one of the swept wings with and without 
a fuselage.^ Steady and unsteady pressure 
distributions are compared with experiments and 
with other calculations made with lifting- 
surface theory. Flutter calculations are 
included for four rectangular wings which are 
essentially identical except for wing thickness. 
The resulting flutter characteristics, are com- 
pared with those obtained by use of lifting- 
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surface theory^ ^ and with experimental 
data. Ip all of the calculations, the num- 
ber and distribution of panels has been varied 
in order to examine convergence of the. results. 


SQUSSA Panel Method 

Description 

The analysis which the SOUSSA program 
implements is based on application of the 
infinite-space Green's function method to the 
fully unsteady linearized velocity-potential 
partial differential equation.^ The result is 
an integral expression for the velocity poten- 
tial at any point in the flow at any time in 
terms of the value of the potential and its nor- 
mal derivative over the sur1;ace of the body and 
its wake. If the field point lies on the sur- 
face of the body, the expression becomes an 
integro-differential delay equation for the 
potential on the surface of the body. Computa- 
tion of the integral over the surface of the 
body is accomplished by surface paneling with 
twisted quadrilateral (hyperboloidal ) panels 
which maintain continuity of the surface, 
although discontinuities in surface slope are 
introduced, The result is a set of 
differential-delay equations in time. The time 
integration can be accomplished directly in the 
time domain by finite-difference procedures.How- 
ever, the Pl.l version of the program3»4 

employs a Laplace-transform solution which 
yields frequency-domain (complex frequency) 
aerodynamics in the form of a matrix equation 
relating the velocity potential (usually 
unknown) to the normalwash distribution (usually 
known from specified shape, motion, and deforma- 
tion of the body). Preinultiplying the matrix 
equation by a matrix relating surface pressures 
or generalized forces to the potential and post- 
multiplying by a matrix relating the normalwash 
to the generalized coordinates permits direct 
calculation of surface pressures or generalized 
forces (weighted integrals of pressure). 

The surface boundary condition implicit in 
the SOUSSA formulation is equivalent to the 
usual no-penetration condition and is automati- 
cally satisifed by the representation obtained 
from the Green's theorem. As a consequence, the 
perturbation potential is zero inside the body 
so that disturbances do not propagate into the 
interior. The far-field boundary condition is 
automatically satisifed by the source and dou- 
blet singularities distributed over the body and 
wake surfaces. 

In the SOUSSA Pl.l code, the velocity 
potential is taken to be constant over each of 
the quadrilateral panels (zeroth-order panels). 
Although higher-order panels have been formu- 
lated for both subsonic and supersonic speeds, 
they have not yet been implemented in the pro- 
gram. The wake in the Pl.l program is assumed 
to have zero thickness and to extend downstream 
approximately in the freestram dirction. It is 
not required to be flat, but its shape remains 
"frozen" during the calculations. 

The SOUSSA Pl.l program structure's^ is 
modular to facilitate inclusion of new capabil- 
ities or improved algorithms and to provide 


restart capability. Conservation of computer 
memory is emphasized to permit application to 
complicated shapes, such as complete aircraft, 
requiring large numbers of panels. Efficient 
computations are possible for multiple frequen- 
cies and/or multiple sets of vibration or defor- 
mation modes because the aerodynamic influence- 
coefficient integrals are independent of both 
mode shape and frequency and because the 
elements of the influence matrix depend on fre- 
quency in a very simple way. 

The Pl.l code employs the data base and 
data-handling utilities of the SPAR finite- 
element structural analysis program.l4 These 
were incorporated because SOUSSA Pl.l was ori- 
ginally intended for the calculation of steady- 
state structural loads and unsteady aerodynamics 
for flutter and gust-response calculation in 
multidisciplinary structural-optimization com- 
putations employing the SPAR structural analy- 
sis. The SPAR components, however, are unnec- 
essary for stand-alone use. More efficient 
methods for stand-alone operation are available. 

SOUSSA Pl.l does not have a built-in geome- 
try preprocessor because it is considered pre- 
ferable for the user to have the flexibility of 
choosing a geometry processor that is appro- 
priate for his needs. 


Program Improvements 

Subsequent to the completion of the initial 
form of the SOUSSA Pl.l code and the publication 
of references 3 and 4, several significant 
improvements have been incorporated and others 
are known to be needed for any future version of 
the program. Among the latter, higher-order 
panels and elimination of SPAR components have 
already been mentioned. Another is transposi- 
tion and revision of the solution algorithm to 
substantially reduce input/output (I/O) opera- 
tions, In the application of steady- or 
unsteady-flow panel methods, such as PAN 
AIR15,16 and SOUSSA, to large problems (hun- 
dreds of panels), the major part of the cost is 
I/O related. Preliminary considerations, based 
on operations count, indicate that reduction of 
I/O operations by transposition, along with the 
other changes mentioned, could reduce the cost 
of large SOUSSA calculations by nearly an order 
of magnitude relative to that for the current 
version of the Pl.l program. Improvements 
already made have reduced the cost of 
calculations by 35 to 50 percent relative to 
that for the initial version. 

Some program modifications already or 
currently being accomplished are described in 
the following subsections. 

Out-of-Core Solver - The computers usually 
used to implement SOUSSA, such as CDC CYBER 170 
series machines, have a usable central memory of 
about 120,000 (decimal) words. When allowance 
is made for the essential code to solve the 
SOUSSA system of algebraic equations relating 
the values of the potential at the panel centers 
to the matrix containing the normalwash vectors, 
the largest complex matrix that can be stored in 
memory is about 250x250. This means, for 
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example, that a wing with unsymmetric loading on 
the upper and lower surfaces could not be panel- 
ed finer than 11x11 panels on each surface, and 
complicated wing/body configurations could not 
be solved in core. 

This limitation was overcome by writing an 
out-of-core solver to solve the system of simul- 
taneous algebraic equations by LU decomposition. 
That is, the coefficient matrix, which is dense, 
nonsymmetrical , and complex, is factored into a 
unit upper triangular matrix U and a lower tri- 
angular matrix L. These matrices are stored on 
disk and the simpler triangular systems of equa- 
tions are solved for each vibration mode of the 
problem. The disk files that are used to hold 
the original matrix and the factored matrix, as 
the latter is created, are both serial files so 
all file positioning is done with reads, 
rewinds, and backspaces, rather than with disk 
random accesses. This has two advantages, it 
makes the code very portable, and it makes it 
easy to estimate I/O efficiency by merely count- 
ing rewinds and backspaces, A drawback of this 
method is that it is slightly less efficient 
than disk random access. 

It is very easy to compute the spectral 
norm condition number of the L and U matrices 
while performing an LU factorization. The con- 
dition number of the solution process is the 
product of these condition numbers and is a mea- 
sure of the amount of error induced by the 
method used to solve the system of simultaneous 
equations. For all cases studied, this was not 
a significant source of error. 

The out-of-core solver was developed 
expressly to permit the use of paneling schemes 
that lead to coefficient matrices too large to 
fit in memory. However, it was found that use 
of the out-of-core solver utilizing a small 
workspace reduces the execution cost of SOUSSA 
Pl.l calculations even when only a moderate num- 
ber of panels are used. The reason for this is 
that the typical computer costing algorithm con- 
tains a term proportional to the product of 
central processor time and central memory space, 
and the time in this term is the time of the 
entire SOUSSA run. 

Kutta Condition - One improvement needed 
but not yet incorporated is related to the Kutta 
condition. The usual requirement of pressure 
continuity across the wake is imposed in the 
calculation of the velocity potential. However, 
the use of constant-potential panels in SOUSSA 
Pl.l prevents accurate direct calculation of 
surface pressures near the trailing edge. The 
present program therefore imposes two conditions 
on the variation of pressure in that region from 
which upper- and lower-surface pressures can be 
evaluated near the trailing edge (See equations 
(5-22) to (5-28) of reference 3.): l)The differ- 
ence between upper- and lower-surface pressures 
is required to vanish as the square root of the 
distance from the trailing edge— a condition 
consistent with lifting-surface theory. 2)The 
average of upper- and lower-surface pressures at 
the trailing edge is taken to be the same as the 
average at the centers of the upper- and lower- 
surface panels adjacent to the trailing edge. 

The first of these conditions is reasonable as 


long as the airfoil tailangle is not large. As 
will be shown later in this paper, however, the 
second condition is not a good approximation and 
should be replaced by a better extrapolation or 
preferably by special trai ling-edge panels. 

This condition was never thought to be an accu- 
rate characterization of pressure behavior near 
the trailing edge, but it was incorporated into 
this first computer implementation as a matter of 
expediency. 


Analytic Make Integration - The SOUSSA Pl.l 
program, as originally coded, used a paneling 
scheme for the wake similar to that for the 
body. It was soon found that for the range of 
reduced frequencies of practical interest, the 
wake would have to be paneled for as many as 
five chord lengths downstream. Also it was 
found that large variations in panel size de- 
grades the numerical accuracy of the current 
version of SOUSSA, because it uses low-order 
elements. This means that the foremost wake 
panels cannot be larger than the panels at the 
trailing edge of the lifting surface. As a con- 
sequence it was often necessary to use more 
panels for the wake than for the oscillating 
lifting surface. Although the computation 
required for the wake panels is simple, the I/O 
service required places the same demand on the 
computer for the wake as it does for the lifting 
surfaces. This is a major part of the execution 
cost of SOUSSA Pl.l, This problem was eased by 
paneling each wake strip with a single panel, 
extending to infinity, and integrating the 
interaction function analytically in the stream 
direction. In effect, this replaced a large 
number of low-order finite elements on each wake 
strip by a single high-order element on which 
disturbances propagate downstream with free- 
stream velocity as a consequence of the require- 
ment that the pressure difference across the 
wake vanish. The streamwise integral over the 
wake is then the same incomplete modified Struve 
function that occurs as part of the unsteady 
kernel of the familiar potential -flow downwash 
integral equation^^ and can be evaluated by 
methods developed for that purpose. The parti- 
cular method used for SOUSSA was the 24-term 
exponential approximation described in reference 
18. In addition to integrating the wake inter- 
action function in the stream direction, it also 
has to be integrated in the cross-stream direc- 
tion. For the paneled wake this integration was 
performed in closed form. However, after per- 
forming the streamwise integration analytically, 
the cross-stream integrand is not integrable in 
closed form, so a Legendre-Gauss quadrature 
formula is used. By varying the quadrature 
order it was found that a very high order was 
needed for convergence. This is due to a singu- 
larity in the complex plane that is near the 
interval of integration for wing panels near the 
trailing edge. 

The strongest part of this near singularity 
occurs in the zeroth-order (steady-flow) term. 
When the cross-stream wake integrand is expanded 
in powers of the frequency, all the terms that 
do not contain the logarithm of the frequency 
can be integrated in closed form. However, only 
the zeroth and first-power terms were integrated 
in closed form. The Legendre-Gauss quadrature 
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was used to integrate only the difference be- 
tween the cross-stream integrand and the zeroth 
and first-power terms of its expansion. This 
can be integrated accurately with a very low 
order quadrature. In SOUSSA Pl.l a four-point 
cross-span quadrature is used over each wake 
strip. Even though this technique of performing 
the spanwise integration is required only for 
wing panels near the trailing edge, the integra- 
tions are performed the same for all wake inte- 
grals to simplify the program. 

Replacing a paneled wake by an analytic 
wake in this manner does not save much CPU 
(central processor unit) time. However, there 
is a very significant saving in I/O cost, memory 
cost, and CPU cost for I/O service. A typical 
SOUSSA Pl.l run with the analytical wake costs 
about half as much as the same run with an ade- 
quately paneled wake. 

There are times when the analytical wake 
cannot be used, for example, if there is another 
lifting surface in the wake. For this reason 
the capability of using the paneled wake was 
left in the program as an option. All of the 
present calculations, however, employed the ana- 
lytical wake integration. 

Evaluation of generalized aerodynamic 
forces - An^i^her improvement of the SOUSSA PI .1 
code that is being developed is more accurate 
calculation of generalized aerodynamic forces. 
After the pressures have been calculated at the 
panel centers, the generalized aerodynamic 
forces (weighted integrals of pressure) are cur- 
rently evaluated by multiplying each value of 
pressure by the associated panel area and by the 
value of the weight function (generalized-force 
mode) at the panel center and summing these pro- 
ducts for all the panels. This process amounts 
to rectangular integration. As will be shown 
later in this paper, it can yield generalized 
aerodynamic forces that converge more slowly 
(with respect to number of panels) than most of 
the component pressures used in their evalua- 
tion. Clearly, a better integration scheme is 
needed. A subprogram to perform the integration 
by Gaussian quadrature (weighted or unweighted 
as required) has been written and is being de- 
bugged but has not been incorporated into SOUSSA 
Pl.l. 

An alternate revision that is being consid- 
ered is integration by parts to permit evalua- 
tion of generalized aerodynamic forces directly 
from the potential. 


Applications and Results 
Aerodynamic Calculations 

The SOUSSA Pl.l code calculates pressures 
at the centers of the surface panels. Since 
these points usually will not be located at 
spanwise stations where experimental data are 
available for comparison, it is necessary to 
interpolate the calculated pressure to theappro- 
priate stations. In this study, the cubic 
spline was used for the spanwise interpolation. 

In the following discussion statements of 
SOUSSA paneling arrays refer to the number of 


panels on one-quarter of the wing surface, e.g., 
the upper right-hand side. A 14-by-22 array, 
for example, indicates 14 panels from leading 
edge to trailing edge and 22 panels from root to 
tip. The left/right symmetry option in SOUSSA 
was used throughout. Since all of the configu- 
rations in this study are vertically symmetric 
(that is, symmetric with respect to the x-y 
plane), the vertical symmetry option of SOUSSA 
was used for calculations at zero angle of 
attack. Thus for a 14-by-22 array of panels, 
the number of unknowns (and hence the order of 
the matrix problem to be solved) is 308 for zero 
angle of attack and 616 for nonzero angle of 
attack. Note that the wake does not introduce 
additional unknowns into the problem. 

Comparison of Pressure Distributions from 
SOUSSA Pl.l and RHOIV - An initial step in eval- 
iTating the^ accuracy of results from the SOUSSA 
Pl.l program is to compare distributions of 
pressure for thin, isolated wings as obtained 
from SOUSSA with values generated by a state-of- 
the-art lifting-surface method such as 
RHOIV^*^. Both steady and unsteady pressures 
for this purpose have been generated for a rec- 
tangular wing of aspect ratio 2.0, and a 45- 
degree swept wing generated by shearing back the 
rectangular one. Both wings have a one-percent 
thick biconvex airfoil. The motion is unit- 
amplitude pitch about an axis through the root- 
trailing-edge point. The Mach number is 0.9, 
and the two reduced frequencies are k = 0.0 and 
0.3. These conditions and the rectangular-wing 
planform were chosen to provide continuity with 
some panel ing-convergence results presented in 
figure 29 and 30 of reference 10. 

For the SOUSSA calculations the half span 
of each wing was divided into a lO-by-10 array 
of panels (Figs. 1 and 2) spaced uniformly 
chordwise, and spaced spanwise in equal incre- 
ments of an angular coordinate defined as the 
inverse cosine of fraction of semi span measured 
from the wing root. This latter distribution is 
called a "cosine distribution." 

The program RHOIV, which implements the 
subsonic kernel -function lifting-surface theory, 
employs a series of lifting-pressure functions 
that are continuous over the wing except at edge 
and hingeline discontinuities, and effects a 
solution by downwash collocation. In the pre- 
sent calculations, 12 collocation points were 
used chordwise at each of 10 span stations. The 
resulting continuous pressure distributions can 
be evaluated everywhere except at discontinu- 
ities. The RHOIV calculations do not account 
for the one-percent wing thickness. 

For the rectangular wing, figures 1(a) and 
(b) show the chordwise distribution of lifting 
pressure coefficient ACp per radian of pitch 
amplitude at the 57.5 percent span station for 
the two frequencies. In figure 1(a) the 
pressure for k=0.0 agree very well except at the 
center of the foremost panel . Such a 
discrepancy is not unexpected where a small num- 
ber of low-order panels (i.e., constant velocity 
potential on each panel) are used in regions 
where potential and pressure are varying 
rapidly. In figure 1(b) for k=0.3, the agree- 
ment of the real -part pressure is about the same 
as for k=0.0. For the imaginary-part pressures. 
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the SOUSSA results are systematically slightly 
larger than the RHOIV results except at the most 
forward panel. 

For the 45-degree-swept wing, the same 10- 
by-10 paneling was retained along with the pitch 
axis through the root trailing edge, and 
pressure distributions are shown in figures 2(a) 
and (b) for the same 57.5 percent span station. 
As can be seen, the differences between the 
SOUSSA and RHOIV pressures are somewhat larger 
than for the unswept wing. 

The observed differences between SOUSSA and 
RHOIV pressures are attributable to two factors: 
(t) wing thickness effects in SOUSSA but not in 
RHOIV; and (ii) SOUSSA results near the leading 
edge are probably not fully conveged with 
respect to the number of panels used. Thickness 
effects on lifting pressure are probably very 
small here because the thickness ratio is only 
one percent. With regard to paneling conver- 
gence, figure 3 of reference 19 shows relevant 
steady-state pressure for a two-dimensional air- 
foil calculated by a zeroth-order panel method 
which for these conditions is comparable to 
SOUSSA Pl.l. The results show that eight panels 
chordwise were sufficient to converge pressures 
over most of the airfoil surface but gave 
pressures that were higher than converged values 
near the leading edge. Those results are con- 
sistent with the SOUSSA/RHOIV comparisons in 
figures 1(a) and 2(a). Paneling convergence for 
SOUSSA results is illustrated and discussed 
later in this paper. (See also Ref. 10). 

Rectangular Wing - Steady and unsteady 
pressure distributions have been calculated by 
SOUSSA Pl.l for an aspect-ratio 3 rectangular 
wing with a five-percent-thick circular-arc air- 
foil for which experimental pressure data are 
available in reference 5. Figure 3 shows the 
wing and the three stations at which pressures 
were measured. These stations are at approxi- 
mately 50, 70, and 90 percent of semi span. Two 
paneling arrays were used in the calculations, 

10 by 10 and 14 by 22. In both arrays the 
panels were uniformly distributed in both direc- 
tions. The steady-state results are compared 
with measured pressures for two angles of attack 
at M = 0.7 in figure 4. As the figures show, 
the SOUSSA results for the two paneling arrays 
are essentially the same except for small 
differences near the wing edges and tip where 
the larger array is expected to resolve the 
pressures more accurately. 

For zero angle of attack (Fig. 4(a)), the 
calculated variations of pressure are in quali- 
tative agreement with experiment, but calculated 
pressure levels are consistently lower than 
experiment. Some of this difference may be 
caused by an effective increase in the model 
thickness due to the boundary layer. The agree- 
ment between calculated and measured pressures 
is, in fact, slightly better at M = 0.7 for 
which the Reynolds number is 4.0 x 10^ and the 
boundary layer is assumed to be thinner than at 
M = 0-24 1 for which the Reynolds number is 
2.2 X 10° (results not shown). The experi- 
mental data, however, show no evidence of flow 
separation. Even though the pressures calcu- 
lated by SOUSSA Pl.l are not in very close 


agreement with experiment, they are in very good 
agreement, over the interior portion of the 
planform, with those obtained from LTRAN3, a 
finite-difference smal 1 -perturbation potential - 
flow code^O (comparison not shown). 

At five degrees angle of attack figure 4(b) 
shows SOUSSA pressures to be in reasonable 
agreement with experimental values except near 
the sharp leading edge where the broad region of 
high suction in the measured upper-surface 
pressures shows evidence of a separation bubble. 
Also, just ahead of the trailing edge the calcu- 
lated pressures (especially those obtained with 
the larger number of panels) level off, that is, 
deviate from the general downward trend, as a 
result of the inadequate implementation of the 
trail ing-edge condition discussed previously. 

The unsteady pressures in figure 5 were 
measured and calculated for the wing oscillating 
in its first bending mode. A lO-by-10 uniform 
array of panels was used in the SOUSSA calcula- 
tions because the steady-state results (Fig. 4) 
showed that to be sufficient to obtain reason- 
ably converged pressures. The resulting surface 
pressures (Fig. 5) are in good agreement with 
the scattered experimental values and with 
values calculated by the kernel -function 
lifting-surface theory. 

Clipped-Tip Delta Wing - The clipped-tip 
delta wing shown in figure 6 is the one used in 
the experimental investigation of reference 6. 
All of the data presented in reference 6 are for 
Mach number 0.9; however, all of the data for 
this wing shown herein are for Mach number 0.4. 
Both sets of data were acquired during the same 
experimental program. 

Steady-state surface pressures calculated 
by SOUSSA Pl.l for the three span stations shown 
in figure 6 are compared with the experimental 
data in figure 6 for a=0° and in figure 7 for 
a = 2°. The paneling arrays used were 8 by 8, 

10 by 10, 12 by 12, and 14 by 14, spaced uni- 
formly both chordwise and spanwise. The curves 
in figure 6 and 7 show that the pressures 
obtained with the 8-by-8 array are essentially 
converged over most of the planform and are in 
good agreement with experiment (for which the 
Reynolds number is 9.0 x 10^). However, the 
pressures near the wing edges converge more 
slowly. Convergence near the edges is improved 
by use of finer paneling locally (e.g., a cosine 
distribution, results for which are not shown). 
The cosine distribution concentrates smaller 
panels near the edges where pressures usually 
have the largest gradients. 

The calculated pressures in figures 6 and 7 
show that immediately ahead of the trailing edge 
the curves tend to level off as they did for the 
rectangular wing. This deviation, as mentioned 
previously, is caused by the condition, still in 
SOUSSA Pl.l, that the average of upper-? and 
lower-surface pressures at the trailing edge is 
taken to be the value calculated at the center 
of the upper- and lower-surface panels adjacent 
to the trailing edge. Also note that in 
figure 7 the upper-surface and lower-surface 
pressures do not close at the trailing edge. The 
failure to close is not a result of the SOUSSA 
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calculation but is caused by the use of cubic- 
spline curves to connect pressure-values obtain- 
ed from SOUSSA Pl.l. 

In addition to the pressures just dis- 
cussed, the convergence of the integrated 
pressures (specifically, the lift-curve slope) 
has also been examined. Reference 4 examined 
the convergence of integrated aerodynamic forces 
as a function of the inverse of N, the number of 
panels in the chordwise or spanwise direction. 
The variation was found to be essentially 
linear, giving confidence to linear extrapola- 
tion toward estimated converged values for N^. 
Figure 8 shows a similar plot of the lift-curve 
slope for the clipped-delta wing as a function 
of the inverse square root of Np, the total 
number of panels in the array, which is an 
equivalent but more general parameter. The 
values shown here also lie along an approximate- 
ly straight line, and the extrapolated estimate 
of the converged value for Np-^o* is indicated. 

The rate of convergence is independent of angle 
of attack, and the value of lift-curve slope 
calculated with the 14-by-14 array (the largest 
number of panels used for this wing) is about 
1.75 percent below the estimated converged 
value. Inasmuch as integration is a smoothing 
process, this is not a particularly impressive 
rate of convergence. In fact, it appears to be 
slower than the convergence of the pressure over 
much of the wing. As was mentioned previously, 
the integration performed in SOUSSA Pl.l to 
generate generalized aerodynamic forces (such as 
lift) from the surface pressures is equivalent 
to rectangular integration. A better integra- 
tion scheme is obviously needed. 

Unsteady pressures have also been measured 
and calculated (with a uniform lO-by-10 array of 
panels) for the wing of figure 6 oscillating in 
pitch about the axis shown. The resulting 
pressures per unit amplitude of oscillation are 
compared in figures 9 and 10 for two frequen- 
cies. The overall agreement is good. The real 
component of the pressure is underpredicted a 
bit over the forward portion of the wing and 
overpredicted over the aft portion where some 
boundary-layer thickening probably reduces Cp 
magnitudes even though the Reynolds number for 
these tests was 9 x 10^ based on average 
chord. At both frequencies there is some evi- 
dence of leading-edge vortex separation at the 
two outboard stations (Figs. 9 and 10(b) and 
(c)) which is, of course, not represented in the 
calculations. 

Swept Tapered Wing With and Without Fuse- 
lage - Steady-state surface pressures have been 
calculated by SOUSSA Pl.l for a wing of aspect 
ratio 4, taper ratio 0.6, quarter-chord sweep 
angle 45^, and NACA 65A006 airfoil. Calcula- 
tions have been made for the wing alone with 16 
panels spanwise distributed uniformly and 22 
panels chordwise also spaced uniformly except 
for subdivision near the leading edge. Calcula- 
tions have also been made for the wing 
with a fuselage (Fig. 11) with a lO-by-10 array 
of panels distributed uniformly on the wing (200 
panels on the half span) and 6 panels circum- 
ferentially around the half fuselage with 34 


panels distributed somewhat irregularly along 
its length, as shown in the figure, for a total 
of 404 panels on the half model. This configu- 
ration approximates the wing-body model of 
reference 7 for which measured pressures are 
available for comparison. The panel configura- 
tion in figure 11 accurately represents the 
wind-tunnel model except that the paneled fuse- 
lage is cylindrical along the length of the wing 
root, and the tunnel sting is not represented. 
The cylindrical fuselage section was a restric- 
tion imposed by the crude geometry preprocessor 
that was used to set up these calculations. 

Pressures for the wing alone and for the 
wing of the wing-body combination are compared 
with experimental values for the wing-body ^ in 
figure 12 for Mach number 0.6 and angle of 
attack 40 . The results for the wing of the 
wing-body combination are in good agreement with 
the experimental values, and those for the wing 
alone are in reasonable agreement with experi- 
ment outboard, as would be expected. However, 
the sharp changes and irregular variation of 
pressure for the wing alone over the subdividied 
panels near the leading edge illustrate one of 
the shortcomings of low-order panel methods, 
that is, sensitivity of the results to paneling 
geometry. Another illustration of this problem 
is seen in the calculated distributions of 
pressure along the fuselage (Fig. 13). In this 
figure pressures are shown along meridian B, 
which is 450 off the vertical, and along 
meridian C which is 75^ off the vertical and 
just above the wing surface (See figure 11). 

The designations B and C are taken from refer- 
ence 7. At zero angle of attack. Figure 13(a), 
the pressures vary smoothly along both meridians 
and agree well with the experimental data, 
except toward the aft end of the fuselage where 
the tunnel sting has not been represented in the 
present calculations. At 4^ angle of attack, 
figure 13(b), the comparison is still good along 
meridian B, but on meridian C a large spurious 
local fluctuation in pressure is calculated ad- 
jacent to the wing- root leading edge and trail- 
ing edge. Since panel lengths on the fuselage 
vary substantially in these regions, the calcu- 
lations were repeated with the two rings of 
panel nodes marked "r" in figure 11 deleted in 
order to see if elimination of the short panels 
adjacent to the leading and trailing edges would 
help to smooth the fluctuations in pressure. 

The results are indicated by the diamond symbols 
in figure 13. Figure 13(b), in particular, 
shows no improvement in the pressure fluctua- 
tions. 

From the beginning of the SOUSSA develop- 
ment it was recognized that introduction of 
higher-order panels (at least linear source and 
quadratic doublet distributions) would be 
required to alleviate the kind of problem shown 
in figure 13(b). 

Flutter Calculations for Rectangular Wings 

Flutter analysis - Reference 12 reports an 
experimental flutter investigation of a series 
of aspect-ratio-5.0 rectangular wings with bi- 
convex (circular arc) airfoils and five thick- 
ness ratios, ranging from about 1.4 to 10 
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percent. Flutter calculations have been made 
for four of these wings with aerodynamics from 
SOUSSA Pl.l and from the FAST kernel -function 
lifting-surface program of reference 11. These 
Galerkin modal analyses were made with natural 
modes of a uniform cantilever beam. Three-mode 
analyses used first and second bending and first 
torsion modes; five-mode analyses added third 
bending and second torsion modes. The three- 
mode and five-mode flutter results are very 
close together, thereby indicating convergence 
with respect to the number of modes used. 

Convergence Analysis - As a preliminary to 
the flutter calculations for comparison with the 
experimental data of reference 12, a study was 
made of the convergence of the SOUSSA Pl.l aero- 
dynamic forces with respect to the fineness of 
the paneling. The configuration studied is the 
rectangular wing of reference 12 with four- 
percent-thick airfoil. The Mach number is 
0.756, and the reduced frequency is 0.114 which 
is close to the reduced frequency at flutter. 

Values of real and imaginary parts of the 
generalized aerodynamic forces are plotted on 
the left-hand side of figure 14 as functions of 
1//Np for the first bending mode (modal index 
1) and the first torsion mode (index 2). The 
calculated values shown by the circle symbols 
are for 8-by-8, lO-by-10, 14-by-14, and 18-by-18 
arrays of panels spaced uniformly both chordwise 
and spanwise. The discontinuities in the ordi- 
nate scale should be noted. The solid curves 
are drawn through the data, and the dashed lines 
are linear extrapolations from the last calcu- 
lated points toward estimated converged values 
for Np^®. For comparison, the square symbols 
on the ordinate axis Indicate values obtained 
from the FAST subsonic-kernel lifting-surface 
analysis of reference 11. The latter are not 
necessarily the values to which the SOUSSA Pl.l 
results are expected to converge because SOUSSA 
Pl.l includes airfoil thickness effects, while 
the lifting-surface analysis does not. These 
convergence results may be compared with those 
of reference 10 which show similar trends. 

On the right-hand side of the figure are 
associated values of the flutter-frequency ratio 
w/wct and the flutter-speed index 
(Again, note the discontinuity of the ordinate 
scale). As on the left, the circles and squares 
indicate SOUSSA and FAST lifting-surface aero- 
dynamics, respectively. The numbers 3 and 5 
indicate three-mode and five-mode flutter calcu- 
lations respectively. The notations 3x3 and 
6x8 indicate the arrangement of downwash col- 
location points used in the FAST calculations. 
The two arrangements produced flutter speeds 
within about one half percent of each other. A 
further check with 8 x 10 collocation points 
gave virtually identical results to those of the 
6x8, and therefore the results are judged to 
be converged with respect to the number of col- 
location points. 

An examination of the real and imaginary 
parts of A22> the weighted twisting moment due 
to first-torsion-mode motion, shows two counter- 
acting effects as the number of panels 
Increases. The real part of A 22 is positive 


and increasing, indicating an increasing diver- 
gent moment that tends to decrease the flutter 
speed. In contrast, the imaginary part of A 22 
is negative and becoming more negative. This 
increase in damping of the first-torsion mode 
tends to increase the flutter speed. The net 
effect is less than one percent decrease of the 
flutter-speed index as paneling increased from 
lO-by-10 to 18-by-18. Although this change is 
small, the indicated trend would bring SOUSSA 
Pl.l results even closer to the lifting-surface 
values if the number of panels were increased 
further. 

Comparison with Experiment - In view of the 
results of the convergence analysis, the SOUSSA 
Pl.l calculations made for comparison with the 
experimental flutter data of reference 12 
employed a uniform 18-by-18 array of panels and 
five vibration modes. The corresponding FAST 
kernel -function calculations also employed five 
vibration modes, and the downwash collocation 
points were arrayed six chordwise by eight 
spanwise at the program-default locations. 

The flutter speeds calculated by the two 
methods for the four-percent-thick wing (Fig. 

15) are close together, and both closely follow 
the experimental trend with Mach number. How- 
ever, they are 10 to 11 percent unconservative 
at the lower Mach numbers up to as much as 14 
percent unconservative at Mach numbers near 
0.9. The flutter frequencies calculated by both 
methods, on the other hand, are very close to 
the experimental values. 

Figure 16 shows measured^^ and calculated 
flutter-speed-index values for four different 
wing thicknesses — 1.4, 4, 6, and 10 percent. 
The two curves for the four-percent-thick wing 
are repeated from figure 15. The solid circles 
for 1.4-percent thickness are slightly below the 
solid squares, and the 6- and 10-percent thick 
results are slightly above, thereby showing a 
small but monotonic effect of thickness. The 
FAST results, of course, contain no aerodynamic 
effects of airfoil thickness. Consequently, the 
differences between FAST results for the . 
different wing thicknesses are caused entirely 
by differences in model mass, stiffness, and 
structural -damping properties. Hence the cal- 
culated aerodynamic effects of thickness at a 
given Mach number are perhaps best assessed by 
comparing the SOUSSA calculations with FAST 
results for the same wing thickness. The FAST- 
SOUSSA differences here range from about one 
percent for the 1.4-percent-thick wing to about 
five percent for the 10-percent-thick wing. 

These differences are small, and the trend 
indicated cannot be confirmed by the measured 
data because of the experimental scatter. 

The new kernel -function results shown in 
figures 15 and 16 do not agree with the kernel - 
function-derived results shown on figure 5 of 
reference 12 which showed close agreement with 
the experiments. The current FAST program is 
believed to be numerically superior to the 
kernel -function program of 1959, and to yield a 
more accurate lifting-surface results. Some of 
the improvements in FAST relative to its pre- 
decessors are discussed in reference 18. 
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A possible explanation of the experimental 
flutter speeds being lower than the calculated 
values is as follows: For thin wings with sharp 

leading edges it has been observed experimen- 
tally that, even at small (nonzero) angles of 
attack (static or dynamic), a localized leading 
edge flow-separation "bubble” occurs. (See also 
the discussion of figure 4(b) above). The 
bubble increases the divergent twisting moment 
above that which would be indicated by attached- 
flow theory. Such an increased moment acts to 
lower the flutter speed. Consequently, it is 
possible that the measured flutter speeds for 
the rectangular wings of figures 15 and 16 and 
reference 12 are lower than those that would be 
calculated by any attached-flow theory. 


Concluding Remarks 

Although the SOUSSA method has the ability 
to handle bodies having arbitrary shapes, 
motions, and deformations, the present study has 
focused on applications to some simple wings in 
steady and unsteady motion, including flutter, 
in order to assess the validity, accuracy, and 
usefulness of the current SOUSSA Pl.l program by 
comparison of results with those of lifting- 
surface theory and existing experimental data. 
SOUSSA Pl.l results are found to be very close 
to those obtained from steady and unsteady 
lifting-surface theory and to agree satisfacto- 
rily with experiment for conditions where agree- 
ment should be expected. 

In addition, the results and experience 
with the program serve to highlight and quantify 
needed improvements which had already been 
recognized. These include use of higher-order 
panels (first-order source, secondorder 
doublet), improved implementation of Kutta con- 
dition, elimination of the SPAR data base and 
data handling utilities, and (based on opera- 
tions count) transposition and revision of the 
solution algorithm. For moderate to large pro- 
blems (a few hundred to many hundreds of panels) 
these changes should reduce the cost of a con- 
verged solution by nearly an order of magnitude. 
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(a) k = 0.0 (b) k = 0.3 

Fig. 1 Calculated lifting-pressure 1 concluded, 

distributions for aspect-ratio-2 
rectangular wing pitching about 
trailing edge. M = 0.9, y = 0.575. 
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(a) k = 0.0 

Fig, 2 Calculated lifting-pressure 

distributions for aspect-ratio-2, 
45-degree-swept wing pitching about root 
trailing edge, M = 0.9, y = 0.575. 


Fig, 3 Rectangular wing model of reference 5. 
All dimensions are in inches. 
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(b) k = 0.3 
Fig. 2 Concluded. 
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(a) a = QO. 

Fig, 4 Steady-state pressure distributions on 
aspect-ratio-3 rectangular wing at M = 0,70. 
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Fig. 5 Unsteady pressures on aspect-ratio-3 
rectangular wing at a = 0°. M = 0.2- 




Fig. 6 Steady-state surface pressure on 

clipped-tip delta wing. M = 0.40, 

cx = oo. 



(c) y = 0.85 
Fig. 7 Concluded. 











Fig. 11 SOUSSA paneling for model of 
reference 7. 





(a) y = 0.20 W y = 0.80 

Fig. 12 Calculated and measured surface Fig. 12 Concluded 

pressures on wing of wing-body 
configuration. M = 0.6, a = 4°. 




(b) a = 4° 

Fig. 13 Concluded. 
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Fig. 14 Paneling convergence of generalized 
aerodynamic forces and flutter 
characteristics for four-percent-thick 
rectangular wing of reference 12. M = 
0.756. 



Fig. 15 Flutter characteristics for four- 
percent-thick rectangular wing of 
reference 12. 


Fig. 16 Effect of thickness on flutter-speed- 
index for the rectangular wings of 
reference 12. 
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